Dynamics and stability of Bose-Einstein condensates with attractive 1/r interaction 
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The time-dependent extended Gross-Pitaevskii equation for Bose-Einstein condensates with at- 
tractive 1/r interaction is investigated with both a variational approach and numerically exact 
calculations. We show that these condensates exhibit signatures known from the nonlinear dy- 
namics of autonomous Hamiltonian systems. The two stationary solutions created in a tangent 
bifurcation at a critical value of the scattering length are identified as elliptical and hyperbolical 
fixed points, corresponding to stable and unstable stationary states of the condensate. The stable 
stationary state is surrounded by elliptical islands, corresponding to condensates periodically oscil- 
lating in time, whereas condensate wave functions in the unstable region undergo a collapse within 
finite time. For negative scattering lengths below the tangent bifurcation no stationary solutions 
exist, i.e., the condensate is always unstable and collapses. 

PACS numbers: 03.75.Kk, 34.20.Cf, 04.40.-b, 05.45.-a 



I. INTRODUCTION 

Bose condensates with atomic interactions, in addition 
to the contact interaction, open the possibility to study 
the properties of degenerate quantum gases in which the 
relative strengths of long- and short-range interactions 
can be continuously adjusted by tuning the contact inter- 
action via a Feshbach resonance. In particular, the Bose- 
Einstein condensation of neutral atoms with electromag- 
netically induced attractive 1 jr interaction has been pro- 
posed. According to O'Dell ct al. 1] "monopolar" quan- 
tum gases could be realized by a combination of 6 appro- 
priately arranged "triads" of intense off-resonant laser 
beams. In that arrangement, the 1/r 3 interactions of the 
retarded dipole-dipole interaction of neutral atoms in the 
presence of intense electromagnetic radiation are aver- 
aged out in the near-zone limit, while the weaker 1/r in- 
teraction is retained. An outstanding feature of this type 
of long-range interaction is that for attractive contact in- 
teraction stable Bose-Einstein condensates are predicted 
that are self-bound (without an additional trap). Despite 
the fact that a self-bound condensate with attractive 1/r 
interaction has not been realized in an experiment so far, 
the physical parameters necessary for creating it exper- 
imentally were discussed in detail by Giovanazzi et al. 
2], who also obtained a theoretical prediction for the 
number of atoms in the self-bound condensate. Further- 
more, the existence of stable monopole and quadrupole 
oscillations around the stationary ground state have been 
demonstrated and analyzed with analytical and numeri- 
cal calculations 

The realization of self-bound condensates with attrac- 
tive 1/r interaction in experiments remains a challeng- 
ing task. However, experimental realization has been 
achieved for a similar system with a long-range inter- 
action, viz. dipolar Bose-Einstein condensates. Here, 
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the long-range dipole-dipole interaction appears next to 
the short-range contact term. The system has attracted 
much attention in recent years [1, [E S S Hi and the 
achievement of Bose-Einstein condensation in a gas of 
chromium atoms with a large dipole moment, has 
opened the way to promising experiments on dipolar 
quantum gases [To| . For example, the experimental ob- 
servation of the dipolar collapse of a quantum gas has 
recently been reported by Koch et al. [TTj], which sets 
in when the contact interaction is reduced below some 
critical value. 

Bosc-Einstcin condensates with attractive 1 jr interac- 
tion exhibit a similar collapse scenario. Collapse of the 
self-bound condensates again only sets in below some 
critical value of the scattering length. It was recently 
shown [12j that these critical values in fact correspond 
to bifurcation points where two solutions of the Gross- 
Pitaevskii equation disappear: one the true ground state, 
the other a collectively excited state. It was also demon- 
strated [l3[ that at the bifurcation the two stationary so- 
lutions exhibit the typical structure known from studies 
of exceptional points (3, ITU [T^ . [TtI [T"8| in open quantum 
systems described by non-Hermitean Hamiltonians. At 
the exceptional points both the energies and the corre- 
sponding wave functions are identical, a situation which 
is forbidden for bound states of the linear Schrodinger 
equation (which have to be orthogonal) but is possible 
here because of the nonlincarity of the Gross-Pitaevskii 
equation. 

Using a complex continuation of the Gross-Pitaevskii 
equation the existence of complex eigenstates at (real) 
negative scattering len gth s below the bifurcation point 
has also been revealed [13|]. The physical interpretation 
of these states is the collapse (or explosion) of the con- 
densate with a decay rate given by the imaginary part of 
the complex eigenvalues of the chemical potential. 

It is the purpose of this paper to analyze self-bound 
spherically symmetric Bose condensates with attractive 
1/r interaction in the vicinity of the bifurcation points 
from the point of view of nonlinear dynamics, and to in- 
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vestigate the time evolution of arbitrary condensate wave 
functions. We do this by solving the time-dependent 
Gross-Pitaevskii equation both by means of a variational 
approach with time-dependent Gaussian wave packets 
and by exact numerical computations using the split- 
operator method. We will show that of the two station- 
ary solutions created in a tangent bifurcation one is dy- 
namically stable and the other unstable, corresponding 
to elliptic and hyperbolic fixed points, respectively. The 
stable state is surrounded by solutions periodically os- 
cillating in time, whereas wave functions in the unstable 
region undergo a collapse within finite time. Below the 
tangent bifurcation no stationary solutions exist, i.e., the 
condensate is always unstable and collapsing. 

The special appeal of investigations of the properties 
of spherically symmetric self-trapped Bose condensates 
with attractive l/r interaction lies in the fact that the 
extremization of the variational mean-field energy can be 
carried out fully analytically. The reason is that in ab- 
sence of a trap the extremization conditions for the mean- 
field energy become quadratic equations, while with a 
trap potential the equations become at least polynomials 
of order 5, caused by the combination of the trap and 
contact interaction terms, regardless of whether or not a 
long-range interaction is present, and of its type. There- 
fore generic properties of Bose condensates can be eluci- 
dated in a very simple and transparent way, and later be 
checked in numerical calculations and in more complex 
systems such as dipolar Bose condensates. 

The evolution of a wave function is determined by 
the extended Gross-Pitaevskii equation for self-trapped 
Bose-Einstein condensates with attractive l/r interaction 
without external trap potential, which reads 



A + 8irNa\ip(r,t)\' 



7rmr,t), (1) 



where the natural "atomic" units introduced in Ref. [12| 
were used. Lengths are measured in units of a "Bohr 
radius" a u — h 2 /(mu), energies in units of a "Rydberg 
energy" E u — h 2 /(2ma 2 ), and times in units of t u — 
h/E u , where u determines the strength of the atom-atom 
interaction [![, and m is the mass of one boson. The 
number of bosons is N . 

The paper is organized as follows. In Sec. [IT] the ba- 
sic equations resulting from a time-dependent variational 
principle applied to the Gross-Pitaevskii equation are 
presented, where the solution of the nonlinear partial dif- 
ferential equation is reduced to a set of ordinary first or- 
der differential equations. The linear stability of the fixed 
points is investigated with variational as well as numer- 
ically exact calculations in Sec. IIIIl In Sec. IIVI the dy- 
namics of the condensate obtained from variational and 
from numerically exact computations is investigated and 
the results are compared. Conclusions are drawn in Sec. 
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II. TIME-DEPENDENT VARIATIONAL 
PRINCIPLE 

We ap ply the time-dependent variational principle 
(TDVP) [ll [H to solving the time-dependent Gross- 
Pitaevskii (GP) equation. Exploiting the scaling prop- 
erties presented in Ref. [l2j and introducing the scaled 
variables 

(f, a, t, i>) = (Nr, N 2 a, N 2 t, N~ 3/2 iP) , (2) 
we bring the system in the form 



i-J(r,t) = H${r,t) 
at 



-A f + Vc + V u 



where the potentials 

V c = 8Tra$(r 1 t)\ 2 , 

./ \r-r'\ 



<P(r,t), (3) 

(4a) 
(4b) 



depend on the coordinates and the wave function, i.e., H 
is a nonlinear operator. The scaling reveals that the sys- 
tem has only one external parameter, viz., a — N 2 a/a u 
[jjj]. If not stated otherwise, we use the scaled variables 
throughout the rest of the paper and omit the tilde in 
what follows. 

An approximate solution ip(t) on a given manifold in 
Hilbert space is obtained by minimizing the quantity 
I = \\i<f>(t) — Hip(t)\\ 2 with respect to </> only, and then 
setting tj) = (j). This means that the time-dependent vari- 
ational principle does not depend on properties of the 
operator H, i.e., the method can be applied to linear and 
nonlinear operators in the same way. Let the wave func- 
tion be parametrized by a set of complex time-dependent 
functions z(t) = {zi(i) , . . . , z n (t)} , i.e., ip(t) = ip(z(t)). 
The TDVP yields a set of ordinary differential equations 
for the motion of z(t) 



Kz 



-ih 



(5) 



with the matrix K and the vector h of the n-dimensional 
system §5§ given by 



K 



dip 
dz 



dz 



h 



dtfj 
Oz 



H 



1> 



(6) 



In time-independent calculations for condensates with 
1 /r interaction Gaussian wave functions have been used 
as an ansatz for the two solutions created in the tangent 
bifurcation 0, G3]. To apply the TDVP to the Gross- 
Pitaevskii equation UJ we choose as a test function a 
radially symmetric Gaussian wave packet 



^( r , t) = e ^ Ar '+^ 



(7) 



with the time-dependent variational parameters z(t) = 
{A(t),j(t)} = {A r (t) + iAi(t),~{ r (t) + *7i(t)}. A similar 
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time-dependent Gaussian trial function has been applied 
in studies of the dynamics of the Gross-Pitaevskii equa- 
tion without 1/r interaction [2l[ . Evaluation of K and h 
in iJSJ) and © yields the ordinary differential equations 



.4 



-4A 2 - L -V 2 , 



(8a) 
(8b) 



where vq and V2 are given as solutions of the two linear 
equations 

(V#H + \W*\il))V 2 = (MVc + V u \iP) , (9a) 

(Mr 2 Wv Q + = (Mr 2 (V c + V u )\i/>) , (9b) 

with 



(^|r 2 ^) 
(V|r 4 |^) 

WW) 
(^|r 2 F c |V) 



^ 3/2 f -27 lA -3/2 


(10a) 


3- 3/ %- 27lA -5/ 2 

8V2 1 ' 


(10b) 


157r3/2 e—^-/ 2 

32V2 


(10c) 


^ 5 / 2 ae- 4 ^4- 3/2 , 


(lOd) 


7T 5/2 e -4 7i ,-5/2 

2 


(lOe) 


37r5/2 ae"^A- 5/2 


(lOf) 


16 6 ^ ' 


(10g) 



Inserting vq and V2 in (j5aj) and (|8b[) leads to the differ- 
ential equations 



A = -4A 2 + 2v / 2 7 re~ 271 f aA t - - 

6 



7 = 6«^4 



Tre- 2 ^ 

2-y/2^i 



(5 - 14aA, 



(11a) 
(lib) 



The imaginary part of Eq. (jllb[) can be integrated ana- 
lytically, 



7i(*) 



4 7T 



(12) 



and guarantees the normalization condition ||^|| = 1. 

The final form of the equations of motion in real num- 
ber representation is given by 



A r = -4(4! - Af) + -j=AT [aA, - 
Ai = 8A r Ai , 



j r = -6Ai + -^VAi (5 - UaAi) , 

\/7T 



(13a) 
(13b) 
(13c) 



which are solved under the initial conditions A r (0) = 
7 r (0) = 0, Ai(0) > for an initially real valued Gaussian 
wave packet. 
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FIG. 1: (Color online) Chemical potentials of the ground state 
(e+) and the nodeless excited state (e_) in the variational 
calculation. They emerge in a tangent bifurcation at a critical 
value of the scaled scattering length a CI — — 37r/8 = —1.1780 
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III. LINEAR STABILITY OF THE 
BIFURCATING STATES 

In this section we study the stability of the two station- 
ary eigenstates that are born in the tangent bifurcation. 
We do this first for the variational solutions, and then 
demonstrate how the essential features carry over to the 
analysis of the numerically exact solutions. 



A. Stability of the variational solutions 

In the time-dependent variational approach, the two 
stationary states of the Gross-Pitaevskii equation ([3]) ap- 
pear as the time-independent solutions (fixed points) of 
the equations of motion (|13[) . Requiring Ai = and 
A r = immediately leads to 








1 ± 




(14a) 
(14b) 



The vanishing of the real part of A implies that the state 
indeed is a stationary Gaussian. The scaled chemical 
potentials e are given by the negative time derivative 
— 7 r of the phase of the wave function in Eq. (|13c[) 



_(var) 



-It = 



5±4,/l + ff 



9tt 



1±i/1 



8« 
3tt 



(15) 



The chemical potentials of the two solutions JT5J) are 
drawn in Fig. [TJ The tangent bifurcation behavior of 
the chemical potential at the critical scattering length 
a cr = — 37r/8 = —1.1780 is clearly visible. The branch 
with the lower chemical potential in Fig. [T] turns out 
to have higher mean-field energy than the other branch. 
Therefore the plus sign refers to the ground state, and the 
minus sign to the collectively excited state. The results 
obtained here via the fixed points of the time-dependent 
equations of motion for the variational parameters fully 
agree with the results derived in Refs. [l|, applying a 
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([TB"]) are calculated with the usual ansatz SA r j(t) — 

6A^e x( For the stationary ground state we find the 
two eigenvalues 



^(var) 



±- 



16i 



1 + 



9tt 



8a 
3tt 



(17) 



-1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 



which are always imaginary for a > —3tt/8 (i.e., above 
the bifurcation point) and correspond to the monopolc 
excitation mode considered in Ref. They describe an 
elliptic fixed point, which is stable. Furthermore, we can 
also give the eigenvalues of the collectively excited state, 



FIG. 2: (Color online) Eigenvalues of the linearization of the 
extended Gross-Pitaevskii system Q for the two stationary 
solutions emerging at the tangent bifurcation. Both the vari- 
ational (var) and the numerically exact (num) solutions are 
shown. The two eigenvalues A+ of the stationary ground state 
(cf. Ref. 0|) are purely imaginary demonstrating that the 
state is an elliptic fixed point. In contrast to this finding, the 
two eigenvalues A_ of the nodeless excited state are purely 
real. One is positive, the other is negative. The state is an 
hyperbolic fixed point. Vanishing real or imaginary parts are 
not shown. 



time- independent variational approach to extremize the 
mean-field energy, and which were used to compare with 
numerically exact solutions of the stationary extended 
Gross-Pitaevskii equation in Refs. [l2l [l3j. 

The linearization of the equations of motion (CCi 
around the stationary solutions A r , Aj are given by 



6A r ± = ±- 



1 + f 2 - 



5A i± = - 



9tt 
32 



1 + fs. ± 1 



9tt 



1+ fa ±1 



~5A r ± . 



(16a) 
(16b) 



The eigenmodes with the eigenvalues A ( - var ) of Eqs. 



. (var) 



±- 



16 



3tt 



"9tt 



1 + fi-l 



(18) 



which are always real for negative scattering lengths 
a > — 37r/8, and one eigenvalue is positive. Hence, the 
eigenstate is unstable and belongs to a hyperbolic fixed 
point. 



B. Stability of the numerically exact solutions 

Numerically exact stationary eigenstates of the Gross- 
Pitaevskii equation JT]) which are the counterpart of the 
variational solutions described above have been calcu- 
lated in Refs. [H, EH- Here, we investigate the stability 
of these states and compare the results with the eigenval- 
ues (fl7|) and (fT8|) . For the stability analysis, we extend 
the procedure applied by Huepe et al. [22j for condensates 
without 1 jr interaction in a harmonic trap to the scaled 
integro-differential equation ([3]). The extended Gross- 
Pitaevskii equation is linearized with the Frechet deriva- 
tive, once the wave function has been decomposed in its 
real ip R (r,t) and imaginary , / (r,i) part. The derivative 
is evaluated at the real- valued stationary states ijj±( r ) 
found in Ref. [12], which leads to two coupled equations, 
viz. 



9 S^ R = -A-e|8m^(r) 2 -2 f dr' t±^-L j 5^(r,t) , (19a) 



dt \ J \r — r 

|^ =(-A-e + 24W±(r) 2 - 2 / dV P^r) S^ R (r, t) + 4V3±(r) / ^ r ^*( r W*(r', *) (1Qb) 
at \ J r - r'\ J \r - r'\ 



For the calculation of the eigenmodes we take the per- (5y/(r,i) = <5?/>o(r)e A " um Note that Si/i R (r) and St/j^r) 
turbations in the form Stp R (r,t) — 6ip R (r)e x( and 
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are not always real- valued. As solutions of the linearized 
equations (|19[) they can become complex wave functions. 

The numerical solution of the linearized wave equa- 
tions is done with the method described in Ref. [l3| . Us- 
ing the "linearized potential" 

W .,^4HW5, ,20, 

J \r-r'\ 



XS^(r) = -S4(r)" - 

XS^ (r) = -8^(r)" - 

U 1 (r)" = --U: i(r)'- 
r 

The wave functions 5ip R (r) and 6ipo(r) are only deter- 
mined up to a complex (normalization) constant, which 
they have in common. Therefore, only two real variables 
are required to set the initial value (5ip R (0), SipQ (0)) = 
(cos a, sin a e 1 ^) for the integration. The real angles 
a and /3, the complex eigenvalue A, and the complex 
initial value = Ui(0) are the parameters which 

have to be searched to find the solutions of the system 
([2Tjl . The initial conditions for the first derivatives are 
(Sip R )' = (Si/j 1 )' = U[ = 0. The correct solutions are 
found when we obtain square integrable wave functions 
and the value of C/ 1 initially assumed is identical to the 
integral 

/>oo 

[/i(0) = 167r/ dr' r'^ ± (r')S^(r') , (22) 
Jo 

which can be calculated once the wave function obtained 
numerically has been determined. 

As the Gross-Pitaevskii equation is a nonlinear differ- 
ential equation the correct normalization of the numer- 
ical wave functions cannot be obtained by a multiplica- 
tion with a normalization constant but is possible by ex- 
ploiting scaling properties of equations (f2~Tj) . In the form 
presented here, the proper normalization and scaling of 
the extended Gross-Pitaevskii system is achieved by the 
scaling factor 

/>oo 

l/v= \\4'\\ 2 =4tt / dr\^{r)\ 2 r 2 (23) 
J o 

and the transformation of all of the following values (cf. 

My- 

W,r,e,t,a,U)-> (v 2 i(),-,v 2 e,^,-^,v 2 U) . (24) 

As can be seen from the system ([2"TT) . the differential 
equations are invariant if in addition to the transforma- 
tions (|24p the eigenvalue is scaled by A — > v 2 \. Thus, 



we can transform (fT9| into a system of three second or- 
der ordinary differential equations, which can be solved 
together with the stationary extended Gross-Pitaevskii 
equation. In order to be in a position to compare with 
the results of the variational approach we assume radi- 
ally symmetric wave functions. The ordinary differential 
equations, which have to be solved, then read 



(21a) 
(21b) 
(21c) 
I 

it is possible to solve the system of differential equations 
(|2ip with the unsealed and non-normalized stationary so- 
lutions i/j± (r) of the nonlinear extended Gross-Pitaevskii 
equation computed numerically if a subsequent scaling of 
the eigenvalues is performed. 

When the solution of the linearized system is carried 
out for the stationary ground state, a pair Ai = — A2 
of purely imaginary eigenvalues are found. These are 
plotted as a function of the scaled scattering length a in 
Fig. O The same calculation for the nodeless station- 
ary state yields a pair of purely real eigenvalues. The 
qualitative agreement between the eigenvalues calculated 
numerically and the results obtained analytically in the 
variational approach is very good, but there are signif- 
icant quantitative differences. Such behavior is known 
from previous studies of the system [12, EH . Similar to 
the case discussed in Ref. [22| , there exist neutral modes 
with the eigenvalue A^ num - ) = 0, and additional imagi- 
nary eigenvalues are found for both stationary solutions. 
Altogether, the numerically exact results confirm the sta- 
bility analysis performed using the analytical eigenvalues 
of the variational approach. 



IV. DYNAMICS OF THE CONDENSATE 

In this section we investigate the time evolution of ar- 
bitrary wave functions. We do this again first by the 
variational approach and then by exact numerical calcu- 
lations. The variational calculations will reveal the dif- 
ferent types of dynamics in the vicinity of the elliptic and 
the hyperbolic fixed points, viz. oscillatory or collapsing 
solutions. Furthermore, the collapse of the condensate 
for a < a CI can be followed as a function of time. The 
numerically exact approach will confirm these findings 
and, in addition, exhibit a larger variety of qualitatively 
different dynamical behavior of the condensate. 



% fl (r)' + (247mV; ± (r) 2 - U(rj) 5^(r) ~ I7i(r)^±(r) , 
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FIG. 3: (Color online) Phase portraits of the dynamics ob- 
tained from time-dependent trajectories of the complex func- 
tion A(t). (a) For a = — 1 the two stationary solutions ap- 
pear as fixed points, (b) They coalesce for a = —1.18. (c) For 
a — —1.3, below the bifurcation point, there are no stationary 
solutions. 



A. Variational approach 



We solve the two ordinary differential equations (113a|) 
and (|13bj) for various scaled scattering lengths a above 
and below the critical value a cr = — 3ir/8 = —1.1780. 
Phase portraits of the dynamics can be obtained by plot- 
ting the imaginary part as a function of the real part of 
the time-dependent trajectories A(t) = A r (t) + iA^if). 
The phase portraits of trajectories with different ini- 



FIG. 4: (Color online) Periodically oscillating condensate for 
the scattering length a = — 1 and the initial condition Ai(0) = 
0.3. 



tial conditions are shown in Fig. [3J for three values of 
the scaled scattering length, one above the bifurcation 
point (a = —1, Fig. [3(a)), one near the bifurcation point 
(a — —1.18, Fig. [3(b)), and one below the bifurcation 
point (a = —1.3, Fig. [3(c)). For a > a CI the elliptic and 
the hyperbolic fixed points are clearly recognizable in Fig. 
[3(a), they correspond to the stationary eigenstates. The 
two fixed points coalesce at the critical scattering length 
a cr (see Fig. [3(b)), and disappear for a < a cr (Fig. [3(c)) 
implying that no longer stationary eigenstates exist. 

For the physical interpretation of the phase portraits 
it is useful to note that the width \J (r 2 )(t) of the con- 
densate is related to Ai (t) via 



(r 2 )(t) 



(Mr 2 m 3 



(25) 



Thus, in the stable region surrounding the elliptic fixed 
point the width of the condensate oscillates periodically. 
This is illustrated in Fig. 2] for a condensate with scat- 
tering length a — — 1 and initial condition Ai(0) = 0.3. 

In the unstable regions Ai(t) increases to infinity, 
which means the collapse of the condensate. In Fig. [5(a) 
the width \J (r 2 ) is shown for a condensate with scaled 
scattering length a = — 1 and initial condition close to 
the hyperbolic fixed point, given by Eq. Ipajh pbl) at 
A r = 0, Ai = 0.3787. The width first stays approxi- 
mately constant, as is to be expected in the vicinity of 
the unstable stationary state, however, as soon as the 
decrease of the width becomes noticeable, the complete 
collapse to zero width occurs within finite time. This be- 
havior demonstrates the existence of a collapse induced 
by the attractive atom-atom interactions. In a realis- 
tic experimental situation further mechanisms, which go 
beyond the scope of this paper, have to be taken into ac- 
count. In particular, the contraction of the wave function 
amplifies density-dependent inelastic collisions which re- 
sult in a loss of particles [23j, [24j and change the time 
evolution. Note that the calculations presented here as- 
sume a constant scaled scattering length (cf. Eq. ([2])). Of 
course, for initial conditions close to the hyperbolic fixed 
point the evolution of the collapse depends sensitively on 
the initial deviation from the fixed point, i.e., the closer 
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(a) 



(b) 




FIG. 5: (Color online) Collapse of the unstable eigenstate. (a) 
Scaled scattering length a = — 1 and initial condition close to 
the hyperbolic fixed point at A r = 0, A% = 0.3787. The 
collapse after finite time is clearly visible in the inset, but 
depends sensitively on the initial conditions (see text), (b) 
Collapse of the non-stationary state at a = —1.3 with initial 
condition Ai(0) = ^ + ^ = 0.10416. 



the initial conditions approach the hyperbolic fixed point 
the longer it takes before the collapse sets in. 

At scattering lengths a < a cr the square root in Eq. 
(I14bp becomes imaginary which means that no fixed 
points and thus no stationary states can exist. The con- 
densate collapses for all initial conditions A(0) of the 
wave function. This is illustrated for a = —1.3 in Fig. 
EJb) with the initial condition of the "least" unstable 
non-stationary state, Ai(0) = ^ + = 0.10416 (cf. 
Eq. (|14b[) ). Contrary to the situation shown in Fig. O^a) 
the decrease of the condensate width starts immediately 
without any plateau. Ignoring a loss of particles we find 
that the width vanishes after the finite time T c = 9.2522. 

The variational approach with complex parameters 
A(t) and jit) in Eq. ([7]) ensures that the mean-field en- 
ergy 



E 



3 (A 2 + A 2 r ) 2^fA~{2aAi - 1) 
Ai 



(26) 



is conserved. Thus, the phase portraits in Fig. [3] can be 
calculated alternatively as equipotential lines E = const 
instead of integrating the equations of motion (|13[) . 

In fact, the equations of motion obtained from the 
time-dependent variational principle can be brought into 



Hamiltonian form if the variational parameters A r , Ai are 
replaced with two other dynamical quantities, of which 
one is assigned to be the momentum and the other the 
coordinate variable. Such adequate canonical variables 
are [H[ 



(27a) 
(27b) 



and in this set the mean-field energy reads 

J? u( \ t^t/ 2,9, 3V3a v/3 
E = H{q,p) = 1 + V = p ' 




4q 2 ^/nq 



(28) 



with the decomposition into a "kinetic" part T depend- 
ing on the "momentum" p and a "potential" part V de- 
pending only on the "coordinate" q. Note that q has the 
physical meaning of the square root of the radius of the 
condensate, according to equation (|25[) . If the mean-field 
energy (|28|) is identified with a Hamiltonian, it describes 
the one-dimensional motion of a particle in the potential 
V(q) obeying Hamilton's equations: 



2p 



dH 
dp 

dH _ 9 
~~dq ~ 2g3 



3 9a 

7T V 



3 1 

7r q 2 



(29a) 
(29b) 



Of course the backward substitution A r = Ai = 
together with (|29a[) and (|29b[) yields the same equa- 
tions of motion for A r and Ai as obtained from the 
TDVP. Conversely, if the trial wave function ([7]) had been 
parametrized by q, p, 7, the TDVP would have yielded 
their equations of motion (|29a[) and (|29b[) . 

The "potential" part V(q) of the mean-field energy 
(|28p is plotted in Fig. [S] as a function of the "position" 
variable q for different scattering lengths below, at, and 
above the critical scattering length a cr . It agrees with the 
mean- field energy plotted in Ref. [l[, calculated using a 
real-valued spherically symmetric (stationary) Gaussian 
trial wave function. In our approach the "kinetic" term 
p 2 in Eq. (f2"8"|) appears additionally in the mean-field en- 
ergy because of the complex ansatz ([7]). In other words, 
the full mean-field energy of Ref. |l| corresponds to our 
potential part V , in which the condensate moves like a 
classical particle. 

Above a cr the potential possesses two stationary 
points, a stable one at the minimum and an unstable 
one at the maximum. At a cr the bifurcation takes place, 
i.e., the two extrema coincide and there is only a saddle 
point of the potential. For a < a cr the potential has no 
stationary points. For a > a cr the motion of the conden- 
sate is stable as long as the mean-field energy lies below 
the maximum of the peak of the potential and if it is lo- 
cated close to the local minimum on the right-hand side 
of the unstable fixed point. The one-dimensional motion 
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FIG. 6: (Color online) "Potential" part of the mean-field en- 
ergy for different scattering lengths. This potential part of 
the mean field energy (|28[l agrees with the complete mean- 
field energy for stationary states of Ref. [lj] (see text). The 
motion of the Gaussian wave function is interpreted as the 
one-dimensional motion of a classical particle in the potential 
V(q). 
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FIG. 7: (Color online) Phase portrait of the mean-field energy 
formulated in the canonical variables q and p for a = —0.8. 



is periodic between two turning points. If the energy is 
increased above the energy of the unstable fixed point, 
only one turning point remains, and the condensate col- 
lapses when q approaches zero. 

Expanding the potential about the minimum and max- 
imum leads to the frequencies (|17p of oscillations of the 
ground state and the decay rates (fT5)) of the collectively 
excited state. The phase space portrait of the dynamics 
in (p, g)-variables is presented in Fig. [7] for the scattering 
length a = —0.8. The equipotential lines which asymp- 
totically approach the q = axis represent the collapse 
of the condensate. 



B. Exact time-dependent calculations with the 
split-operator method 

The split-operator method [26| assumes the decompo- 
sition of the Hamiltonian H = T + V, and makes use 
of the short-time approximation for the time evolution 
operator 

e - lT( T+V) = e -i(T/2)T e -iTV e -i(T/2)T + 0(r 3) _ qq) 

The kinetic part of the evolution operator is applied to 
the wave function in momentum space, the potential part 
is applied to the wave function in position space rep- 
resentation. The method is numerically especially effi- 
cient when a fast Fourier transform is used for the tran- 
sition from momentum to position space representation 
and backwards. For the nonlinear Gross-Pitaevskii equa- 
tion ([3]) the potential part of the time evolution opera- 
tor needs some further investigation in comparison with 
the linear Schrodinger equation. Although the scatter- 
ing term V c belongs to the nonlinear part of the Gross- 
Pitaevskii equation, it can be treated like a conventional 
potential of a linear Schrodinger equation and presents 
no additional difficulty The 1 /r interaction potential V u 
however requires the additional solution of the integral 
(|4b[) after every time step of integration. This integral is 
computed using the convolution theorem, i.e., 

HVu(r,t)} = H~}HMr,t)\ 2 } , (31) 
r 

where the Fourier transform of 1 jr is performed analyt- 
ically to give 

,1, 4t , s 

H~} = — ■ 32 

r p z 

Altogether two additional fast Fourier transforms are 
necessary per time step and we obtain for the 1/r in- 
teraction term 

1 (\ Z" 00 * poo 

Vj r ,t) = -— dp==f- dr'r'Mr',t)\ 2 smpr' . 
r Jo P Jo 

(33) 

The fast Fourier transforms are performed on equidis- 
tant grids with 1024 or 2048 points. The numerical con- 
vergence is checked by monitoring that the wave function 
vanishes at the r max , p max border of the grid in position 
as well as in momentum space. Due to the widening 
of the wave function in some computations the neces- 
sary size of the grid depends on the desired propagation 
time. A criterion of convergence for the time step of the 
split-operator method is that the mean-held energy is a 
constant of motion. In particular for the computations 
in which the wave function is very close to the hyperbolic 
fixed point the time step must be chosen small for conver- 
gence. Specifically, for the computations in the vicinity 
of the hyperbolic fixed point we use At = 0.0001 whereas 
for computations far away from the unstable stationary 
state a time step of At = 0.01 is sufficient. 
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In the following paragraphs, the dynamics of the con- 
densate is investigated for regions initially close to the 
two stationary states. In particular, we present the evo- 
lution of those initial states which are obtained by de- 
forming the stable and the unstable stationary state by 



ip-n/>-f, r~*r/f 



2/3 



(34) 



= 1.48 at t = 
increasing time and 



with a stretching factor /, i.e., for / = 1 the stable and 
the unstable stationary state, respectively, are obtained. 
This choice of perturbation leaves the norm of the wave 
function unchanged. 

We start with the investigation of the dynamics of wave 
functions close to the unstable stationary state. In Fig. 
Ufa) the square root of the width (|2"5")l of the condensate 
is plotted as a function of time. The evolution presented 
in Fig.JSfa) is computed at the scaled negative scattering 
length a = —0.85 with a deformation factor / = 1.001. 
The wave function itself is plotted in Fig. [5](b) as a func- 
tion of the radial coordinate r at different times. The 
initial wave function is given by the solid red line. The 
condensate stays nearly stationary at the beginning for 
times up to t sa 4, i.e., the dashed green line in Fig. 
[H](b) representing the wave function at t = 4 nearly co- 
incides with the initial state, and the width (Fig. Ufa)) 
has only slightly decreased from \J (r 2 
to yp)" = 1.44 at t = 4.0. With 
the farther away the wave function has moved from the 
stationary state, the shrinking of the width of the wave 
function is accelerated. As already predicted by the vari- 
ational computation this leads to a collapse of the con- 
densate at finite time. Both Figs. [51(a) and [FJb) show 
that the width of the condensate tends to zero in posi- 
tion space when the collapse time T c is approached. Con- 
versely, in momentum space the wave function becomes 
arbitrarily wide close to T c . Thus, the size of the grid 
in momentum space is the numerically limiting factor for 
the propagation with the split-operator method. Choos- 
ing a large grid in momentum space allows for integrating 
arbitrarily close to T c . 

According to the variational results there exist peri- 
odic solutions in the vicinity of the unstable stationary 
state, and, indeed, similar behavior is found in the nu- 
merically exact computations. In Fig. HJa) the width of 
a quasi-periodically oscillating condensate at a = —1.0 
is shown. The associated wave function is plotted in 
Fig. E0(b) for different times. It can be seen that at the 
times where the root-mean-square extension of the con- 
densate goes through a maximum, or minimum, there is 
also a good agreement of the respective wave functions. 
This shows that the wave functions indeed oscillate quasi- 
periodically. The calculation is started with the wave 
function of the unstable stationary state (/ = 1). Be- 
cause of numerical deviations, the solution begins to os- 
cillate for times larger than f w 25. In contrast to the 
variational result, the oscillation of the condensate is not 
strictly periodic here. 

The dynamics of the condensate presented so far does 
not qualitatively differ from what is predicted by the vari- 




FIG. 8: (Color online) (a) Width of a slightly perturbed sta- 
tionary state (a — —0.85, / = 1.001) as a function of time, 
(b) Wave functions of the state for selected times marked by 
circles in (a). The collapse of the condensate is obvious from 
both plots. 



ational calculation. In Fig. fiWah however, the situa- 
tion is encountered where the width increases monoton- 
ically with time. Here we have chosen a deformation of 
/ = 0.99 at the scaled scattering length a = —0.85. For 
short times there is again a plateau where the width is 
nearly constant as in Fig. [9ja). The plateau is shorter 
than that in Fig. OJa) because the initial wave function 
of Fig. fiWb) differs more from the stationary state than 
the initial state associated with Fig. OJa). For times 
t > 65, however, the width increases linearly with time. 
Obviously Fig. [TUTb) shows a broadening of the wave 
function for times up to t « 20. For longer times the main 
peak of the wave function located at the origin does not 
broaden monotonically with time as might be concluded 
from Fig. fTOT a). Indeed, the wave functions at t = 240 
and t = 400 seem to be even sharper than the wave func- 
tion at t = 20 in Fig. fTUT b). This apparent contradiction 
is resolved, if the wave functions are plotted on a broader 
range in position space and on a logarithmic scale, which 
is presented in Fig. [TTJ Here it is clearly visible that the 
wave functions have more than one peak and their width 
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FIG. 9: (Color online) Evolution of the unstable "stationary" 
state (/ = 1) at a = —1.0. (a) Due to numerical deviations 
an oscillation of the width starts for times larger than t ~ 25. 
(b) Corresponding wave functions for the times marked in (a). 



does increase for larger times. However, the amplitude 
of the run-away parts is very small compared to the first 
maximum at the origin. For obtaining the accurate prop- 
agation it is therefore necessary to choose a large grid in 
position space to make the wave function vanish on the 
border, in particular for long propagation times. The 
usual method to prevent the wave function from run- 
ning to the border of the grid and being reflected there, 
namely to introduce an absorbing complex potential, is 
not applicable for the nonlinear Gross-Pitaevskii system, 
in which an absorption of the wave function alters the 
potential. 

In the vicinity of the stable stationary state the con- 
densate exhibits an oscillatory behavior shown in Fig. 
I12f a). The parameter values for this computation are a 
stretching of / = 1.01 of the stable stationary state at 
a = —0.85. The amplitude of the oscillation is very small 
compared to the oscillation in Fig. 0(b) and ranges from 
about yj (r 2 ) = 3.336 to yj (r 2 ) = 3.385. This observation 
is in correspondence with the variational result in Fig. 
inja) where small deviations from the stable fixed point 
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FIG. 10: (Color online) (a) The width of the condensate, 
initially in the vicinity of the unstable fixed point with / = 
0.99 and a = —0.85 increases linearly with time for t > 65. (b) 
Wave functions of the same state for different times marked 
by circles in (a). The wave functions at t — 240 and t — 400 
seem to be sharper than the wave function at t = 20. Note 
that a larger step size on the grid in position space than in 
Figs. [5] and was used. Therefore there is a visible distance 
between r = and the first point on the grid. 



lead to small oscillations around the fixed point along 
the equipotcntial lines whereas small deviations from the 
unstable fixed point may lead to oscillations with a large 
amplitude. However, with increasing distance from the 
stable fixed point the dynamics is no longer oscillatory, 
but, as can be seen in Fig. [12jb), there is also a gradual 
broadening of the condensate, which becomes more pro- 
nounced for larger deviations of the initial wave function 
from the stable stationary state as is shown in Fig.[T2Tb), 
in which the time evolution of the width for / = 1.25 is 
plotted. The scattering length is set to a = —0.85 again. 
Here, the oscillatory motion at the beginning changes to 
a mainly expanding motion with minor modulations. 
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FIG. 11: (Color online) Double logarithmic plot of the wave 
functions presented in Fig. llOf b - ). The long range tail oc- 
curring with increasing propagation times and leading to a 
monotonically increasing width in Fig. HOf a) becomes visible. 



V. CONCLUSION 

We have investigated the time-dependent extended 
Gross-Pitaevskii equation for self-bound Bose-Einstein 
condensates with attractive l/r interaction. The two sta- 
tionary states of the system, which emerge in a tangent 
bifurcation, and were found with time-independent cal- 
culations in Refs. [l], [13], were identified as an elliptic 
and a hyperbolic fixed point of the dynamical system. A 
time-dependent variational approach opened the possibil- 
ity to obtain analytical results for the eigenvalues of the 
linearized equations of motion, which is a feature of the 
self-trapped condensate with l/r interaction. Numeri- 
cally exact calculations confirmed the variational find- 
ings. 

The dynamics of the system was analyzed with the 
time-dependent variational approach and numerically ex- 
act computations with the split-operator method. For 
scaled scattering lengths larger than the critical value an 
oscillatory behavior of the condensate was found in the 
vicinity of the stationary solutions. It was shown that it 
is possible to introduce canonical variables and interpret 
the dynamics of the width of the condensate as the mo- 
tion of a classical particle in a potential. Furthermore, 
there are unstable regions of the phase space, in which a 
collapse of the condensate occurs. For scaled scattering 
lengths below the critical value the condensate is always 
collapsing. The numerically exact method revealed solu- 
tions with a continuously expanding motion of the width. 

The self-bound spherically symmetric Bose-Einstein 
condensates with gravity-like interaction investigated in 
this paper are unique in so far as variational calculations 
to a large extent can be performed analytically. The re- 
sults obtained in this way can serve as a useful guide to 
numerical calculations, and to the exploration of the non- 
linear dynamic properties of more general Bose-Einstein 
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FIG. 12: (Color online) (a) Quasi periodically oscillating 
condensate at a scaled scattering length of a — —0.85 and 
/ = 1.01. (b) An initially quasi periodically oscillating con- 
densate at a — —0.85, / = 1.25 turns into expanding dynam- 
ics after long time. 

systems. The appearance of stable elliptic islands, with 
elliptic fixed points, unstable regions, with hyperbolic 
fixed points, and of separatrices dividing the regions, are 
typical signatures of autonomous Hamiltonian systems. 
Therefore they will also be present if a symmetric trap 
potential is added to the contact and gravity-like inter- 
action. 

In anisotropic condensates the variational approach 
will lead to autonomous Hamiltonian dynamics with 
more than one degree of freedom. Therefore in addi- 
tion to the signatures discussed in this paper signatures 
of chaos can appear. In fact, a recent variational investi- 
gation of Bose-Einstein condensates with a dipole-dipole 
interaction [27( has revealed that, as the mean-field en- 
ergy is increased from its ground-state value, chaotic re- 
gions emerge in phase space which coexist with regular 
islands, corresponding to irregularly fluctuating or quasi- 
periodically oscillating states of the condensates, respec- 
tively. The investigations presented in this paper have 
paved the way for studies of this sort in realistic Bose- 
Einstein condensates with long-range interactions. 
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